clear all
set more off

global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"

cd $data

cap log close
log using $figures\K_file,replace t


************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
************************************************************************************************************************************************************************************

global controls_LV 	female college racex2 doctor* married widowed age  walkra dressa stoopa beda hosp bmi i.occ_interm experience
global controls_V 	$controls_LV vign_f i.vign_id

global controls_L 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.occ_interm i.wave 
global controls_A 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.occ_interm i.wave logbs 
global controls_R 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.occ_interm i.year bsx1-bsx9

global controls_C 	female college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave ///
					i.occ_interm logy lowcash govt_hins priv_hins 		
global controls_0C	$spec_C spin somecov* sp_medexp*	

************************************************************************************************************************************************************************************
************************************************************************************************************************************************************************************


*******BEGIN PROGRAMS TO BE USED IN BOOTSTRAP

***To get vignette parameters
cap program drop vignette_reg
program define vignette_reg
	qui biprobit (vign_severe = $controls_V) (cantwork =$controls_LV) ,iterate(30)
	gen bF_V	=_b[vign_severe:female]
	gen bFV_V	=_b[vign_severe:vign_f]
	gen rho_V_L	=r(table)[1,colsof(r(table))]
	gen cvg_VL= e(converged)
	matrix bV=e(b)
	matrix bV=bV'
end

***ML Program for the joint estimation of {L,A,R} 
cap program drop myll2
program define myll2
args lnf xb1 xb2 xb3 c21 c31 c32
        tempvar sp2 sp3 k1 k2 k3
        quietly {
            gen double `k1' = 2*$ML_y1 - 1
            gen double `k2' = 2*$ML_y2 - 1
            gen double `k3' = 2*$ML_y3 - 1
            tempname cf21 cf31 cf32 cf22 cf33 C1 C2
            su `c21', meanonly
            scalar `cf21' = r(mean)
            su `c31', meanonly
            scalar `cf31' = r(mean)
            su `c32', meanonly
            scalar `cf32' = r(mean)
            scalar `cf22' = sqrt( 1 - `cf21'^2 )
            scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^2 )
            tempname C1 C2
            mat `C1' = (1, 0, 0 \ `cf21', `cf22', 0 \ `cf31', `cf32', `cf33')
            mat `C2' = (1, 0 \ `cf21', `cf22' )
            egen `sp3' = mvnp(`xb1' `xb2' `xb3') if $ML_y1==1, ///
                    chol(`C1') dr(10) prefix(z) signs(`k1' `k2' `k3') adoonly
            egen `sp2' = mvnp(`xb1' `xb2' ) if $ML_y1==0, ///
                    chol(`C2') dr(10) prefix(z) signs(`k1' `k2') adoonly
            replace `lnf'= ln(`sp3') if $ML_y1==1
            replace `lnf'= ln(`sp2') if $ML_y1==0
        }
end

***Program tp predict Xb_L and Xb_A (stat discrim)
cap program drop predict_a1a2
program define predict_a1a2
	#delimit ;
	gen a1_fromV=
	bV[2 ,1]		*college+  
	bV[3 ,1]         *racex2  +
	bV[4 ,1]   *doctor_told_has_hbp +
	bV[5 ,1]   *doctor_told_has_psy +
	bV[6 ,1]   *doctor_told_has_hea +
	bV[7 ,1]   *doctor_told_has_art +
	bV[8 ,1]   *doctor_told_has_dia +
	bV[9 ,1]   *doctor_told_has_lun +
	bV[10,1]   *doctor_told_has_str +
	bV[11,1]   *doctor_told_has_can +
	bV[12,1]         *married +
	bV[13,1]         *widowed +
	bV[14,1]             *age +
	bV[15,1]      *walkra_R_1 +
	bV[16,1]      *dressa_R_1 +
	bV[17,1]      *stoopa_R_1 +
	bV[18,1]        *beda_R_1 +
	bV[19,1]          *hosp_R +
	bV[20,1]          *bmi_R  +
	bV[21,1]     *(occ_interm==1)+
	bV[22,1]     *(occ_interm==2)+
	bV[23,1]     *(occ_interm==3)+
	bV[24,1]     *(occ_interm==4)+
	bV[25,1]     *(occ_interm==5)+
	bV[26,1]     *(occ_interm==6)+
	bV[27,1]     *(occ_interm==7)+
	bV[28,1]     *(occ_interm==8)+
	bV[29,1]      *experience;
	#delimit cr
	#delimit ;
	gen a1_inst_fromV=
	bV[2 ,1]		*college+  
	bV[4 ,1]   *doctor_told_has_hbp +
	bV[5 ,1]   *doctor_told_has_psy +
	bV[6 ,1]   *doctor_told_has_hea +
	bV[7 ,1]   *doctor_told_has_art +
	bV[8 ,1]   *doctor_told_has_dia +
	bV[9 ,1]   *doctor_told_has_lun +
	bV[10,1]   *doctor_told_has_str +
	bV[11,1]   *doctor_told_has_can +
	bV[14,1]             *age +
	bV[15,1]      *walkra_R_1 +
	bV[16,1]      *dressa_R_1 +
	bV[17,1]      *stoopa_R_1 +
	bV[18,1]        *beda_R_1 +
	bV[19,1]          *hosp_R +
	bV[20,1]          *bmi_R  +
	bV[21,1]     *(occ_interm==1)+
	bV[22,1]     *(occ_interm==2)+
	bV[23,1]     *(occ_interm==3)+
	bV[24,1]     *(occ_interm==4)+
	bV[25,1]     *(occ_interm==5)+
	bV[26,1]     *(occ_interm==6)+
	bV[27,1]     *(occ_interm==7)+
	bV[28,1]     *(occ_interm==8)+
	bV[29,1]      *experience;
	#delimit cr
#delimit;
	gen a1=
	_b[L:female]*				      female+
	_b[L:college]*                   college+
	_b[L:age]*                           age+
	_b[L:racex2]*                     racex2+
	_b[L:experience]*             experience+
	_b[L:married]*                   married+
	_b[L:widowed]*                   widowed+
	_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[L:doctor_told_has_psy]* doctor_told_has_psy+
	_b[L:doctor_told_has_hea]* doctor_told_has_hea+
	_b[L:doctor_told_has_art]* doctor_told_has_art+
	_b[L:doctor_told_has_dia]* doctor_told_has_dia+
	_b[L:doctor_told_has_lun]* doctor_told_has_lun+
	_b[L:doctor_told_has_str]* doctor_told_has_str+
	_b[L:doctor_told_has_can]* doctor_told_has_can+
	_b[L:walkra_R_1]*             walkra_R_1+
	_b[L:dressa_R_1]*             dressa_R_1+
	_b[L:stoopa_R_1]*             stoopa_R_1+
	_b[L:beda_R_1]*                 beda_R_1+
	_b[L:bmi_R]*                       bmi_R+
	_b[L:hosp_R]*                     hosp_R+
	_b[L:1b.wave]*                  (wave==1)+
	_b[L:2.wave]*                   (wave==2)+
	_b[L:3.wave]*                   (wave==3)+
	_b[L:4.wave]*                   (wave==4)+
	_b[L:5.wave]*                   (wave==5)+
	_b[L:6.wave]*                   (wave==6)+
	_b[L:7.wave]*                   (wave==7)+
	_b[L:8.wave]*                   (wave==8)+
	_b[L:9.wave]*                   (wave==9)+
	_b[L:10.wave]*                  (wave==10)+
	_b[L:11.wave]*                  (wave==11)+
	_b[L:12.wave]*                  (wave==12)+
	_b[L:13.wave]*                  (wave==13)+
	_b[L:14.wave]*                  (wave==14)+
	_b[L:15.wave]*                  (wave==15)+
	_b[L:1b.occ_interm]*        (occ_interm==1)+
	_b[L:2.occ_interm]*         (occ_interm==2)+
	_b[L:3.occ_interm]*         (occ_interm==3)+
	_b[L:4.occ_interm]*         (occ_interm==4)+
	_b[L:5.occ_interm]*         (occ_interm==5)+
	_b[L:6.occ_interm]*         (occ_interm==6)+
	_b[L:7.occ_interm]*         (occ_interm==7)+
	_b[L:8.occ_interm]*         (occ_interm==8)+
	_b[L:_cons];
	#delimit cr
	replace a1=. if L==.

	#delimit;
	gen a2=
	_b[A:female]*				      female+
	_b[A:college]*                   college+
	_b[A:age]*                           age+
	_b[A:racex2]*                     racex2+
	_b[A:experience]*             experience+
	_b[A:married]*                   married+
	_b[A:widowed]*                   widowed+
	_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[A:doctor_told_has_psy]* doctor_told_has_psy+
	_b[A:doctor_told_has_hea]* doctor_told_has_hea+
	_b[A:doctor_told_has_art]* doctor_told_has_art+
	_b[A:doctor_told_has_dia]* doctor_told_has_dia+
	_b[A:doctor_told_has_lun]* doctor_told_has_lun+
	_b[A:doctor_told_has_str]* doctor_told_has_str+
	_b[A:doctor_told_has_can]* doctor_told_has_can+
	_b[A:walkra_R_1]*             walkra_R_1+
	_b[A:dressa_R_1]*             dressa_R_1+
	_b[A:stoopa_R_1]*             stoopa_R_1+
	_b[A:beda_R_1]*                 beda_R_1+
	_b[A:bmi_R]*                       bmi_R+
	_b[A:hosp_R]*                     hosp_R+
	_b[A:1b.wave]*                  (wave==1)+
	_b[A:2.wave]*                   (wave==2)+
	_b[A:3.wave]*                   (wave==3)+
	_b[A:4.wave]*                   (wave==4)+
	_b[A:5.wave]*                   (wave==5)+
	_b[A:6.wave]*                   (wave==6)+
	_b[A:7.wave]*                   (wave==7)+
	_b[A:8.wave]*                   (wave==8)+
	_b[A:9.wave]*                   (wave==9)+
	_b[A:10.wave]*                  (wave==10)+
	_b[A:11.wave]*                  (wave==11)+
	_b[A:12.wave]*                  (wave==12)+
	_b[A:13.wave]*                  (wave==13)+
	_b[A:14.wave]*                  (wave==14)+
	_b[A:15.wave]*                  (wave==15)+
	_b[A:1b.occ_interm]*        (occ_interm==1)+
	_b[A:2.occ_interm]*         (occ_interm==2)+
	_b[A:3.occ_interm]*         (occ_interm==3)+
	_b[A:4.occ_interm]*         (occ_interm==4)+
	_b[A:5.occ_interm]*         (occ_interm==5)+
	_b[A:6.occ_interm]*         (occ_interm==6)+
	_b[A:7.occ_interm]*         (occ_interm==7)+
	_b[A:8.occ_interm]*         (occ_interm==8)+
	_b[A:logbs]*logbs+
	_b[A:_cons];
	#delimit cr
	replace a2=. if A==.

	#delimit;
	gen a3=
	_b[R:female]*				      female+
	_b[R:college]*                   college+
	_b[R:age]*                           age+
	_b[R:racex2]*                     racex2+
	_b[R:experience]*             experience+
	_b[R:married]*                   married+
	_b[R:widowed]*                   widowed+
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020)+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8)+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9+
	_b[R:_cons];
	#delimit cr
	replace a3=. if R==.
end


***To get Spending/Limitation&Application parameters
cap program drop C_A_reg
program define C_A_reg
	qui heckman logx $controls_C, select(nonzero= $controls_0C ) nolog iterate(30) 
	gen bF_C=_b[logx:female]
	gen sigma_C=r(table)[1,colsof(r(table))-1]
	gen cvg_C= e(converged)
	predict xb,xb
	gen u=logx-xb
	qui biprobit (uncensored_C =$controls_0C) (L =	$controls_L),nolog	iterate(30)
	gen cvg_0L= e(converged)
	predict auc1,xb1
	predict auc2,xb2
	scalar rho_0L=e(rho)			/****Don't use scalar!!!!*/ /*???*/
	gen corr_C_0=normalden(auc1)*(normal((1-rho_0L^2)^(-0.5)*(auc2-rho_0L*auc1)))*binormal(auc1,auc2,rho_0L)^(-1)	
	gen corr_C_L=normalden(auc2)*(normal((1-rho_0L^2)^(-0.5)*(auc1-rho_0L*auc2)))*binormal(auc1,auc2,rho_0L)^(-1)
	qui reg u corr_C_0 corr_C_L if L==1 & uncensored_C==1,nocons
	gen E_uC_L=_b[corr_C_L]
	drop auc1 auc2 corr_C*
	qui biprobit (uncensored_C =$controls_0C) (A =	$controls_A),nolog	iterate(30)
	gen cvg_0A= e(converged)
	predict auc1,xb1
	predict auc2,xb2
	scalar rho_0A=e(rho)
	gen corr_C_0  =normalden(auc1)*(normal((1-rho_0A^2)^(-0.5)*(auc2-rho_0A*auc1)))*binormal(auc1,auc2,rho_0A)^(-1)			
	gen corr_C_App=normalden(auc2)*(normal((1-rho_0A^2)^(-0.5)*(auc1-rho_0A*auc2)))*binormal(auc1,auc2,rho_0A)^(-1)
	qui reg u corr_C_0 corr_C_App if A==1 & uncensored_C==1,nocons
	gen E_uC_A=_b[corr_C_App]
	scalar rho_C_0_TBUB		=_b[corr_C_0]
	scalar rho_C_App_TBUB	=_b[corr_C_App]
	gen Aw=1-R
	qui biprobit (uncensored_C =$controls_0C) (Aw =	$controls_R),nolog iterate(30)
	gen cvg_0Aw_R= e(converged)
	predict auc3,xb1
	predict auc4,xb2
	scalar rho_unc_aw=e(rho)
	egen c2=rmean(auc1 auc3)
	gen  c3=A_App_3ml
	gen  c4=A_Aw_3ml
	scalar r23=rho_0A
	scalar r24=rho_unc_aw 
	scalar r34=rho_aw_ap
	scalar r34_2=(r34-r23*r24)/(sqrt(1-r23^2)*sqrt(1-r24^2))
	scalar r24_3=(r24-r23*r34)/(sqrt(1-r23^2)*sqrt(1-r34^2))
	scalar r23_4=(r23-r24*r34)/(sqrt(1-r24^2)*sqrt(1-r34^2))

	gen cvg_partcorr = (r34_2^2<=1) & (r24_3^2<=1) & (r23_4^2<=1)

	scalar r34_2=max(min(r34_2,1),-1)							/*Constrain the partial corr. coeff. to be b/w -1 and 1*/ 
	scalar r24_3=max(min(r24_3,1),-1)
	scalar r23_4=max(min(r23_4,1),-1)

	matrix Sigma = (1, r23, r24 \ r23, 1, r34 \ r24, r34, 1)

	matrix symeigen X lamda = Sigma

	gen cvg_chol = (lamda[1,1]>0) & (lamda[1,2]>0) & (lamda[1,3]>0)

	forvalues j=1(1)3	{
		matrix lamda[1,`j']=max(lamda[1,`j'],0.000001)			/*Replaces negative eigenvectors with a small positive number to avoid error of n.p.d. matrix*/
	}
	matrix Sigma2=X*diag(lamda)*X'
	mdraws, neq(3) dr(20) prefix(q) random

	matrix C=cholesky(Sigma2)
	egen F3 = mvnp(c2 c3 c4), chol(C) dr(20) prefix(q) adoonly
	gen corr_C_00=normalden(c2)*binormal((1-r23^2)^(-0.5)*(c3-r23*c2),(1-r24^2)^(-0.5)*(c4-r24*c2),r34_2)/F3
	gen corr_C_ap=normalden(c3)*binormal((1-r23^2)^(-0.5)*(c2-r23*c3),(1-r34^2)^(-0.5)*(c4-r34*c3),r24_3)/F3
	gen corr_C_aw=normalden(c4)*binormal((1-r24^2)^(-0.5)*(c2-r24*c4),(1-r34^2)^(-0.5)*(c3-r34*c4),r23_4)/F3
	gen v=u-(rho_C_0_TBUB)*corr_C_00-(rho_C_App_TBUB)*corr_C_ap
	qui reg v corr_C_aw if Aw==1 & uncensored_C==1 & A==1,nocons		
	gen E_uC_R=_b[corr_C_aw]
end


cap program drop bayes_invert
program define bayes_invert
	scalar rho_euA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)								/*We want the corr b/w epsilon and uA, not between uL and uA*/
	gen pDi=binormal(xA,xD,rho_euA)/normal(xA)
	gen m1i=normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,xD,rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,xD,rho_euA)^(-1)	
	gen m0i=-normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,-xD,-rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(1-normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,-xD,-rho_euA)^(-1)	
	scalar sigma=sqrt(1+s_csi^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)
	qui su part1 if female==1 & A==1
	scalar p1_1=r(mean)
	qui su part1 if female==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	qui su part2 if female==1 & A==1
	scalar p2_1=r(mean)
	qui su part2 if female==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	qui su part3 if female==1 & A==1
	scalar part3m=r(mean)
	probit R female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9
	predict PI,pr
	gen psi=invnorm(PI)
	gen rho_FX		=part1+part2
	gen lambda_FX	=part3*female
	su psi 		if female==1 & R!=.
	scalar psi1=r(mean)
	su psi 		if female==0 & R!=.
	scalar psi0=r(mean)
	su rho_FX 		if female==1 & R!=.
	scalar rho1=r(mean)
	su rho_FX 		if female==0 & R!=.
	scalar rho0=r(mean)
	su lambda_FX 	if female==1 & R!=.
	scalar lambda1=r(mean)
end
	
*******END PROGRAMS TO BE USED IN BOOTSTRAP
cd $data
	
*FIRST DRAW OF BOOTSTRAP	
*************************************************************************************************	
*************************************************************************************************	
u vignette_data, clear
set seed 2134569
bsample,cluster(hhidpn)

*&*	
	vignette_reg
*&*

keep bF_V bFV_V rho_V_L cvg*
keep in 1
gen draw=1
sort draw
save V_pars,replace


*************************************************************************************************	
cd $data
u rejection_data, clear
set seed 2134569
bsample,cluster(hhidpn)

// Create Halton draw variables
mdraws, dr(10) neq(3) prefix(z)			/*If dr(N), then in egen `sp3' and egen `sp2' adjust dr(N) accordingly*/

gen A=applied
gen R=nosuccess
gen L=cantwork_notemp

// Get initial values
quietly {
        probit A $controls_A
        mat b1 = e(b)
        mat coleq b1 = A
        probit L $controls_L
        mat b2 = e(b)
        mat coleq b2 = L
        probit R $controls_R
        mat b3 = e(b)
        mat coleq b3 = R
        mat b0 = b1, b2, b3
}

ml model lf myll2 	(A: A=$controls_A) ///
					(L: L=$controls_L) ///
					(R: R=$controls_R) ///
					/c21 /c31 /c32 ///
					, missing

ml init b0
qui ml maximize, iterate(30)

gen cvg_ALR= e(converged)

scalar betaFL=_b[L:female]	
scalar delta =_b[R:female]	/*to be used for pinning down \alpha*/

gen bF_L=_b[L:female]
gen bF_A=_b[A:female] 
gen bF_R=_b[R:female]
	
predict_a1a2
gen  A_App_3ml=a2
gen  A_Aw_3ml=-a3

nlcom (r21: /c21) (r31: /c31) (r32: /c21*/c31 + sqrt( 1 - /c21^2 )*/c32) 

scalar rho_aw_ap=-r(b)[1,2]	  
gen rho_A_L=r(b)[1,1]
gen rho_A_R=r(b)[1,2]
gen rho_L_R=r(b)[1,3]

preserve
	gen xbL=(a1-betaFL*female)-(a1_inst_fromV)				
	keep hhidpn a1 a2 xbL female A R college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R year occ_interm bsx1-bsx9
	save $data\data_to_est_alpha_bootstrap,replace			
restore	  
	  
drop if wave<3 			
duplicates tag hhidpn year logx,gen(check_type12_only)

gen nonzero=x>0
gen uncensored_C=logx!=.

*&*	  
	C_A_reg	  
*&*
	
keep in 1
gen draw=1
drop if draw==.
keep draw rho_* cvg* bF_* E_uC_* sigma_C
sort draw
save R_pars,replace

u R_pars,clear
merge draw using V_pars

	**********
		keep bF_V bFV_V rho_V_L bF_L bF_A rho_A_L rho_A_R rho_L_R bF_C sigma_C E_uC_L E_uC_A E_uC_R
		order bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L  E_uC_R E_uC_A rho_A_R
		mkmat bF_L- rho_A_R,mat(M)
		mat M1=M'
		local endmata end
		mata: mata clear
		cap mata: mata drop eval()
		cap mata: mata drop get_baseline_pars
		cap mata: mata drop get_baseline_pars.eval0()
		mata:
		void eval0(todo, b, v, g, H)
		 {
				  x1 =(b[4]-b[5])/sqrt(1+b[1]^2+b[9]^2)										
				  x2 =(b[4]-b[5]-b[6])/sqrt(1+b[1]^2+b[9]^2+b[2]^2)							
				  x3 =(1+b[1]^2+b[9]^2)/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2+b[2]^2))	
				  x4 =-(b[5]-b[7])/sqrt(1+b[1]^2+b[9]^2)											
				  x5 =(b[4])/sqrt(1+b[1]^2+b[9]^2)											
				  x6 =b[9]^2/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2))					
				  x7 =-1/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[11]^2))								
				  x8 =b[8]*b[4]+b[10]															
				  x9 =sqrt(b[8]^2+b[3]^2)														
				  x10=b[8]/sqrt(1+b[1]^2+b[9]^2)												
				  x11=b[8]/sqrt(1+b[11]^2)
				  x12=b[8]/sqrt(1+b[1]^2+b[9]^2+b[2]^2)
				  x13=-1/(sqrt(1+b[1]^2+b[9]^2+b[2]^2)*sqrt(1+b[11]^2))
				  M1=st_matrix("M1")
				  M2=st_matrix("M2")
				  m1 =M1[1,1]
				  m2 =M1[2,1] 
				  m3 =M1[3,1]
				  m4 =M1[4,1]
				  m5 =M1[5,1] 
				  m6 =M1[6,1]
				  m7 =M1[7,1] 
				  m8 =M1[8,1]
				  m9 =M1[9,1]
				  m10=M1[10,1]
				  m11=M1[11,1]
				  m12=M1[12,1]
				  m13=M1[13,1]
			v=((x1-m1)^2/1)+((x2-m2)^2/1)+((x3-m3)^2/1)+((x4-m4)^2/1)+((x5-m5)^2/1)+((x6-m6)^2/1)+((x7-m7)^2/1)+((x8-m8)^2/1)+((x9-m9)^2/1)+((x10-m10)^2/1)+((x11-m11)^2/1)+((x12-m12)^2/1)+((x13-m13)^2/1)
		 }
		S = optimize_init()
		optimize_init_evaluator(S, &eval0())
		optimize_init_evaluatortype(S, "v0")
		optimize_init_params(S, (1,1,1,1,1,1,1,1,1,1,1))
		optimize_init_which(S, "min")
		optimize_init_which(S, "min")
		optimize_init_conv_ptol(S, 1e-16)
		optimize_init_conv_vtol(S, 1e-16)
		b = optimize(S)
		b'
		st_matrix("beta",b')
		`endmata'

		scalar s_omega=beta[1,1]
		scalar s_v=beta[2,1]
		scalar s_phi=beta[3,1]
		scalar mu=beta[4,1]
		scalar mu=beta[5,1]
		scalar gamma=beta[6,1]
		scalar theta=beta[7,1]
		scalar lamda=beta[8,1]
		scalar s_psi=beta[9,1]
		scalar psi=beta[10,1]
		scalar s_csi=beta[11,1]
		scalar sigma	=sqrt(1+s_csi^2)
	**********
		
u data_to_est_alpha_bootstrap,clear
gen xD=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*female			/*Assumption is that X's enter L*, not L_bar*/
ren a2 xA

*&*
	bayes_invert
*&*

gen alpha1=(delta-part1m-part2m)/part3m
gen alpha1_alt=((psi1-psi0)-(rho1-rho0))/lambda1
	
keep alpha1*
keep in 1
gen draw=1
su
save A_pars,replace
*************************************************************************************************	

*BOOTSTRAP DRAWS FROM 2 TO 208
*************************************************************************************************	
forvalues j=2(1)208	{
	u vignette_data, clear
	set seed 2`j'34569
	bsample,cluster(hhidpn)
		vignette_reg
	keep bF_V bFV_V rho_V_L cvg*
	keep in 1
	gen draw=`j'
	sort draw
	save V_pars_`j',replace	
	u V_pars,clear
	append using V_pars_`j'
	save,replace


	cd $data
	u rejection_data, clear
	set seed 2`j'34569
	bsample,cluster(hhidpn)

	// Create Halton draw variables
	mdraws, dr(10) neq(3) prefix(z)			/*If dr(N), then in egen `sp3' and egen `sp2' adjust dr(N) accordingly*/

	gen A=applied
	gen R=nosuccess
	gen L=cantwork_notemp

	// Get initial values
	quietly {
			probit A $controls_A
			mat b1 = e(b)
			mat coleq b1 = A
			probit L $controls_L
			mat b2 = e(b)
			mat coleq b2 = L
			probit R $controls_R
			mat b3 = e(b)
			mat coleq b3 = R
			mat b0 = b1, b2, b3
	}

	ml model lf myll2 	(A: A=$controls_A) ///
						(L: L=$controls_L) ///
						(R: R=$controls_R) ///
						/c21 /c31 /c32 ///
						, missing

	ml init b0
	qui ml maximize, iterate(30)

	gen cvg_ALR= e(converged)

	scalar betaFL=_b[L:female]	
	scalar delta =_b[R:female]	/*to be used for pinning down \alpha*/

	gen bF_L=_b[L:female]
	gen bF_A=_b[A:female] 
	gen bF_R=_b[R:female]
		
	predict_a1a2
	gen  A_App_3ml=a2
	gen  A_Aw_3ml=-a3
		
	nlcom (r21: /c21) (r31: /c31) ///
		  (r32: /c21*/c31 + sqrt( 1 - /c21^2 )*/c32)

	scalar rho_aw_ap=-r(b)[1,2]	  
	gen rho_A_L=r(b)[1,1]
	gen rho_A_R=r(b)[1,2]
	gen rho_L_R=r(b)[1,3]

	preserve
		gen xbL=(a1-betaFL*female)-(a1_fromV)
		keep hhidpn a1 a2 xbL female A R college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R year occ_interm bsx1-bsx9
		save $data\data_to_est_alpha_bootstrap,replace			
	restore	  
		  
	drop if wave<3 			
	duplicates tag hhidpn year logx,gen(check_type12_only)

	gen nonzero=x>0
	gen uncensored_C=logx!=.
		  
		C_A_reg	  

	keep in 1
	gen draw=`j'
	drop if draw==.
	keep draw rho_* cvg* bF_* E_uC_* sigma_C
	sort draw
	save R_pars_`j',replace

	u R_pars,clear
	append using R_pars_`j'
	save,replace

	u R_pars_`j',clear
	merge draw using V_pars_`j'
	
	*************THIS IS EWMD (with I as weighting matrix)
		keep bF_V bFV_V rho_V_L bF_L bF_A rho_A_L rho_A_R rho_L_R bF_C sigma_C E_uC_L E_uC_A E_uC_R
		order bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L  E_uC_R E_uC_A rho_A_R
		mkmat bF_L- rho_A_R,mat(M)
		mat M1=M'
		local endmata end
		mata: mata clear
		cap mata: mata drop eval()
		cap mata: mata drop get_baseline_pars
		cap mata: mata drop get_baseline_pars.eval0()
		mata:
		void eval0(todo, b, v, g, H)
		 {
				  x1 =(b[4]-b[5])/sqrt(1+b[1]^2+b[9]^2)										
				  x2 =(b[4]-b[5]-b[6])/sqrt(1+b[1]^2+b[9]^2+b[2]^2)							
				  x3 =(1+b[1]^2+b[9]^2)/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2+b[2]^2))	
				  x4 =-(b[5]-b[7])/sqrt(1+b[1]^2+b[9]^2)											
				  x5 =(b[4])/sqrt(1+b[1]^2+b[9]^2)											
				  x6 =b[9]^2/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2))					
				  x7 =-1/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[11]^2))								
				  x8 =b[8]*b[4]+b[10]															
				  x9 =sqrt(b[8]^2+b[3]^2)														
				  x10=b[8]/sqrt(1+b[1]^2+b[9]^2)												
				  x11=b[8]/sqrt(1+b[11]^2)
				  x12=b[8]/sqrt(1+b[1]^2+b[9]^2+b[2]^2)
				  x13=-1/(sqrt(1+b[1]^2+b[9]^2+b[2]^2)*sqrt(1+b[11]^2))
				  M1=st_matrix("M1")
				  M2=st_matrix("M2")
				  m1 =M1[1,1]
				  m2 =M1[2,1] 
				  m3 =M1[3,1]
				  m4 =M1[4,1]
				  m5 =M1[5,1] 
				  m6 =M1[6,1]
				  m7 =M1[7,1] 
				  m8 =M1[8,1]
				  m9 =M1[9,1]
				  m10=M1[10,1]
				  m11=M1[11,1]
				  m12=M1[12,1]
				  m13=M1[13,1]
			v=((x1-m1)^2/1)+((x2-m2)^2/1)+((x3-m3)^2/1)+((x4-m4)^2/1)+((x5-m5)^2/1)+((x6-m6)^2/1)+((x7-m7)^2/1)+((x8-m8)^2/1)+((x9-m9)^2/1)+((x10-m10)^2/1)+((x11-m11)^2/1)+((x12-m12)^2/1)+((x13-m13)^2/1)
		 }
		S = optimize_init()
		optimize_init_evaluator(S, &eval0())
		optimize_init_evaluatortype(S, "v0")
		optimize_init_params(S, (1,1,1,1,1,1,1,1,1,1,1))
		optimize_init_which(S, "min")
		optimize_init_which(S, "min")
		optimize_init_conv_ptol(S, 1e-16)
		optimize_init_conv_vtol(S, 1e-16)
		b = optimize(S)
		b'
		st_matrix("beta",b')
		`endmata'

		scalar s_omega=beta[1,1]
		scalar s_v=beta[2,1]
		scalar s_phi=beta[3,1]
		scalar mu=beta[4,1]
		scalar mu=beta[5,1]
		scalar gamma=beta[6,1]
		scalar theta=beta[7,1]
		scalar lamda=beta[8,1]
		scalar s_psi=beta[9,1]
		scalar psi=beta[10,1]
		scalar s_csi=beta[11,1]
		scalar sigma	=sqrt(1+s_csi^2)
	**********

	erase R_pars_`j'.dta
	erase V_pars_`j'.dta
		
	u data_to_est_alpha_bootstrap,clear
	gen xD=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*female			/*Assumption is that X's enter L*, not L_bar */
	ren a2 xA

		bayes_invert

	gen alpha1=(delta-part1m-part2m)/part3m
	gen alpha1_alt=((psi1-psi0)-(rho1-rho0))/lambda1
	
	keep alpha1*
	keep in 1
	gen draw=`j'
	save A_pars_`j',replace

	u A_pars,clear
	append using A_pars_`j'
	erase A_pars_`j'.dta
	save,replace
	}

*************************************************************************************************	
cd $data

u V_pars,clear
sort draw
save,replace
u R_pars,clear
sort draw
save,replace
u A_pars,clear
sort draw
save,replace

u V_pars,clear
sort draw
merge draw using R_pars
drop if _merge!=3
drop _merge
sort draw
merge draw using A_pars
drop if _merge!=3
drop _merge

egen nonconv=rsum(cvg_VL cvg_ALR cvg_C cvg_0L cvg_0A cvg_0Aw_R cvg_partcorr cvg_chol)
keep if nonconv==8			/*Only keep bootstraps where all regressions converge*/

set seed 110968
gen uu=uniform()
sort uu
keep in 1/200				/*Randomly keep 200 bootstrap replications out of the 208 ran*/

drop uu draw

ren alpha1 alpha_old
ren alpha1_alt alpha1		/*This is the DD method*/

su alpha1
scalar se_alpha1=r(sd)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R alpha1,cov 
matrix V=r(C)

ren E_uC_L rho_C_L
ren E_uC_A rho_C_A
ren E_uC_R rho_C_R
keep rho*
collapse (sd) rho*
save $data\sd_bootstrap,replace

/*This comes from J_file: Needs point estimates from final row of Table 6 (corr. coeff.)*/
u $data\moments_for_R,clear		
keep if namevar=="corr_L_A"|namevar=="corr_L_V"|namevar=="corr_L_R"|namevar=="EuC_L"|namevar=="EuC_Aw"|namevar=="EuC_A"|namevar=="corr_R_A"
replace M1=-M1 if namevar=="EuC_Aw"
replace namevar="EuC_R" if namevar=="EuC_Aw"
mkmat M1
u $data\sd_bootstrap,clear
mkmat  rho_V_L-rho_C_R
matrix define z=rho_A_L\rho_V_L\rho_L_R\rho_C_L\rho_C_R\rho_C_A\rho_A_R
matrix z=M1,z
matrix rownames z = corr_L_A corr_L_V corr_L_R corr_C_L corr_C_R corr_C_A corr_R_A 
matrix colnames z = estimate std_error 
esttab matrix(z) using $figures\table6_2.tex, cells("mean(fmt(%9.3f))") label replace style(tex) noobs nodepvar nonumber


log close
clear


